#!/usr/bin/perl
# 
# Exercise 9.5
# 
# Include a graphic display of the cut sites in the restriction map by printing the 
# sequence and labeling the recognition sites with the enzyme name. Can you make a 
# map that handles multiple restriction enzymes? How can you handle overlapping 
# restriction sites?
#
# Answer to Exercise 9.5 - Part 2

use strict;
use warnings;
use BeginPerlBioinfo;

#
# To handle multiple restriction sites, we'll have to make a few changes.  First,
# we'll allow the user to type in multiple restriction enzyme names, and gather
# them into an array.  Then, we'll have to redo the part of the code that
# calls formatrestrictionmap, so that instead of having one subroutine that does
# everything, we'll have a few subroutines:
#
# initialize_annotation:
#    that just creates a blank annotation string of the same length as the sequence; then
#
# add_annotation:
#    that adds the annotation for a given enzyme and location array
#
# formatrestrictionmap2
#    that prints out the annotation and the sequence.

# Declare and initialize variables
my %rebase_hash = ();

unless(dbmopen(%rebase_hash, 'DBNAME', 0644)) {

	print "Cannot open DBM file DBNAME with mode 0644\n";
}

my @file_data = ();
my $query = '';
my $dna = '';
my $recognition_site = '';
my $regexp = '';
my @locations = ();

# If there is a command-line argument, assume it's DNA
if(@ARGV) {
    $dna = $ARGV[0];

# Otherwise, prompt for a FASTA file to open
}else{

    print "Input a FASTA filename: ";
    my $filename = <STDIN>;
    chomp $filename;

    unless(-e $filename) {
    	print "$filename does not exist!\n";
	exit;
    }
    # Read in the file 
    @file_data = get_file_data($filename);

    # Extract the DNA sequence data from the contents of the FASTA file
    $dna = extract_sequence_from_fasta_data(@file_data);
}

# Get the REBASE data into a hash, from file "bionet.110"
%rebase_hash = parseREBASE('bionet.110');

#
# Prompt user for restriction enzyme names, create restriction map
#
do {
    print "Search for what restriction site or sites (or quit)?: ";
    
    $query = <STDIN>;

    chomp $query;

    # Exit if empty query
    if ( ($query =~ /^\s*$/) or ($query =~ /^\s*quit/i) ) {

        exit;
    }

    # Make an array of the restriction site names
    my @queries = split(' ', $query);

    # Initialize the annotation string
    my $annotation = initialize_annotation($dna);

    foreach my $query (@queries) {

        # Perform the search in the DNA sequence

	# If the requested restriction enzyme $query is not known
        if ( not exists $rebase_hash{$query} ) {
	    print "The restriction enzyme $query is not known!\n";
	    next;
	}

	# Obtain the regular expression of the restriction site
        ($recognition_site, $regexp) = split ( " ", $rebase_hash{$query});

        # Create the restriction map
        @locations = match_positions($regexp, $dna);

	# If the restriction sites do not exist in the sequence
	if( not @locations) {
            print "A restriction site for $query is not in the DNA:\n";
	}
    
        # Report the presence of the restriction site and the locations
	#  of the restriction map to the user
        print "Searching for $query $recognition_site $regexp\n";
        print "A restriction site for $query at locations:\n";
        print join(" ", @locations), "\n";
    
	# Create the annotation of the restriction sites
    	$annotation = add_annotation($annotation, $query, @locations);
    
        print "\n";

    }
    
    # Format the restriction map as sequence and annotation
    my @output = formatrestrictionmap2(50, $dna, $annotation);

    print @output;
    
} until ( $query =~ /quit/ );

exit;

################################################################################
# Subroutines
################################################################################


################################################################################
#
# Subroutine initialize_annotation
#
#  Make a blank string of the same length as the given sequence string 

sub initialize_annotation {

    my($seq) = @_;

    return ' ' x length($seq);
}


################################################################################
#
# Subroutine add_annotation
#
#   Add annotation to an annotation string
#

sub add_annotation {

    my($annotation, $enz, @pos) = @_;

    # Put the labels for the enzyme name at the correct positions in the annotation
    foreach my $location (@pos) {
        substr($annotation, $location-1, length($enz)) = $enz;
    }

    return $annotation;
}

################################################################################
#
# Subroutine formatrestrictionmap2
#
# Interleave a sequence with annotation lines showing the locations of an enzyme
#

sub formatrestrictionmap2 {

    my($line_length, $seq, $annotation) = @_;

    # Make a formatted array of the sequence
    my @formatted_sequence = format_sequence($seq, $line_length);

    # Make a formatted array of the annotation
    my @formatted_annotation = format_sequence($annotation, $line_length);

    # Merge the formatted annotation with the formatted sequence

    my @output = ();

    while(@formatted_annotation) {
        push(@output, shift @formatted_annotation);
        push(@output, shift @formatted_sequence);
	# separate the lines
	push(@output,"\n");
    }

    # Return the merged annotation and sequence
    return @output;
}
    

sub format_sequence {

    my($seq, $line_length) = @_;

    my(@output) = ();

    # Split sequence into lines of $line_length
    for ( my $pos = 0 ; $pos < length($seq) ; $pos += $line_length ) {
        push(@output, substr($seq, $pos, $line_length) . "\n");
    }

    return @output;
}

